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Abstract 

Numerical continuation methods for deterministic dynamical systems have been 
one of the most successful tools in applied dynamical systems theory. Continuation 
techniques have been employed in all branches of the natural sciences as well as in 
engineering to analyze ordinary, partial and delay differential equations. Here we show 
that the deterministic continuation algorithm for equilibrium points can be extended to 
track information about metastable equilibrium points of stochastic differential equa- 
tions (SDEs). We stress that we do not develop a new technical tool but that we 
combine results and methods from probability theory, dynamical systems, numerical 
analysis, optimization and control theory into an algorithm that augments classical 
equilibrium continuation methods. In particular, we use ellipsoids defining regions of 
high concentration of sample paths. It is shown that these ellipsoids and the distances 
between them can be efficiently calculated using iterative methods that take advantage 
of the numerical continuation framework. We apply our method to a bistable neural 
competition model and a classical predator-prey system. Furthermore, we show how 
global assumptions on the flow can be incorporated - if they are available - by relating 
numerical continuation, Kramers' formula and Rayleigh iteration. 

Keywords: Numerical continuation, bifurcation analysis, metastability, stochastic dy- 
namics, covariance, Lyapunov equation, ellipsoids, iterative methods, neural competition, 
predator-prey system, Rayleigh iteration, Kramers' law. 



1 Introduction 

Consider a deterministic dynamical system given by a differential equation 

^ = x' = L(x;ri (1) 

where x represents phase space variables, /x G K is a parameter and L is an operator or a 
map that describes a deterministic equation e.g. an ordinary differential equation (ODE), 
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partial differential equation (PDE) or delay differential equation (DDE). Time-independent 
solutions of ((TJ) with x' = are steady states (or equilibria) x* = x*(/j) with L(x*(fi); //) = 0. 
Given an equilibrium x*(fii), numerical continuation allows us to efficiently compute how it 
changes under parameter variation i.e. to compute x*(/j,2) for small — /^l- In the case of 
an ODE we have a vector field 

L(x; = f(x; n) with / : R n x R -> R n . 

Numerical continuation can be used to compute a curve of equilibrium points 7 = {x = 
which solves the algebraic equations f(x;fi) = 0. Furthermore, one can compute 
so-called test (or bifurcation) functions for each point on this curve that indicate a change 
of stability of the equilibrium point under parameter variations. 

Introductions to numerical continuation can be found in [29J [TTj, H5j H]. There are also 
many software packages available with various standard continuation algorithms and test 
functions such as MatCont [H 123, AUTO [3Ql [65], PyDSTool [22] and DDE-BIFTOOL 
[32] . The literature on the applications of numerical continuation techniques is extremely 
large. For example, it can be used to compute periodic and homoclinic orbits [71], stable 
and unstable invariant manifolds of equilibrium points |64| . slow manifolds [47] and canard 
orbits [26] in fast-slow systems as well as isochrons [85J, just to name a few. Application 
areas range from physics [Ml [80], chemistry [611 [28] and biology [721 1115] to engineering 
[103, 92J. It is even possible to implement continuation methods directly in experiments 
[TUB] . 

Despite this success story, there seems to be very little work to extend continuation ideas 
to stochastic differential equations (SDEs). Current numerical approaches to SDEs mostly 
focus on simulation and forward integration j59j[8T]. Other available methods are set- valued 
techniques [25] to track invariant measures and the direct solution of forward or backward 
Kolmogorov PDEs [110L 1102] . An approach that tries to utilize classical continuation for 
stochastic problems is the moment map formulation [HI [33] where the primary motivation 
seems to arise from equation-free modelling [TTj . 

However, suppose we have already used numerical continuation for a deterministic ODE 
and found stable equilibrium points or more general stable invariant sets. Then it is a nat- 
ural question to ask how small noise influences the stability of these objects. In general, 
we expect a change to metastable invariant sets [51 [15] so that noise-induced transitions be- 
tween different stable sets can occur. In this paper, we show that there is a very natural and 
straightforward extension of equilibrium continuation in the context of SDEs that provides 
local information about metastable equilibrium points. Our approach can be applied during 
a numerical continuation calculation or, slightly less efficiently, as a post-processing tool. 

Remark: We note that the algorithm we develop here is expected to extend to much 
wider classes of problems such as nonstationary solutions [68] as well as SPDEs [91 J and 
SDDEs [S3]. 

The method is based on combining well-known results and numerical techniques from 
different areas of mathematics and computing. A reader interested in getting an overview of 
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our main steps should consider the analytical example presented in Section [2j The general 
development based on minimal local assumptions is presented in Sections IBB We test our ap- 
proach for a planar vector field modelling neuronal competition in Section [71 In this example, 
we focus on the algorithmic performance and show how to integrate the algorithm in stan- 
dard numerical continuation software. In Section [8] we consider the Rosenzweig- Mac Arthur 
predator-prey system and demonstrate that important dynamical systems conclusions and 
direct interpretations for applications can be obtained from our computational framework. 
Further examples of how our algorithm relates to important conclusions regarding the dy- 
namics of a system can be found in [161 02] ■ m Section [9] a special case with a global 
gradient-structure assumption is considered. 

Preliminary Remark 1: All computations have been carried out in MatLab [78J, version 
R2010b on a standard quad-core 2.4 GHz CPU with 4 GB RAM. The numerical continuation 
calculations of deterministic equilibrium points use version 2.5.1. of cLMatCont |44j . 

Preliminary Remark 2: All norms refer to the Euclidean norm so that we simply use 
the notation || - || instead of || ■ 1 1 s _ All vectors are assumed to be column vectors. The 
superscript notation ( ) T will denote the transpose of vectors/matrices and I is going to 
denote an identity matrix of suitable size for the algebraic operation considered. 

2 An Analytical Example 

We start with a well-known analytical example to motivate the type of problems we are 
interested in and to present the basic conceptual ideas for the numerical analysis. Consider 
the following 1-dimensional SDE with additive noise 

dx t = (/ix t — xf) dt + adW t =: f(x t ; (j,)dt + adW t (2) 

where W t is standard Brownian motion [84] . a controls the noise level and /i 6 1 is the 
main bifurcation parameter. Systems of the form (J2]) appear very frequently in applications 
ranging from mean-field and Ising-type models for phase transitions in classical physics 
[EEIED], reaction-rate theory in chemistry [53 [IBS], single neuron modelling [75] and bistable 
ecosystems [50] in biology. The deterministic part of the SDE is a normal form for a pitchfork 
bifurcation [TH H3]. The dynamics of ()2|) is easily understood by writing it as a gradient 
system 

dx t = -VU^{x)dt + adW t with U^x) := -^x 2 + ^x 4 . (3) 

so that the stochastic process x t can be interpreted as a particle moving in a potential 
Un(x). There is always one trivial deterministic equilibrium for fl3]) given by x = x* = 0. 
For \i < the equilibrium x* is globally attracting and corresponds to a unique minimum of 
the potential U^(x). At /i = a pitchfork bifurcation occurs; see Figure [U The equilibrium 
x* is destabilized and becomes a local maximum (saddle point) of the potential and two 
new locally stable equilibria x ± = ±^//I appear for \x > corresponding to minima of 
Ufj,(x). Although we can easily obtain the deterministic equilibrium curves given in Figured] 
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Figure 1: Bifurcation diagram (top) for the deterministic part of (j2])-([3]). The potentials are 
shown as well (bottom). In the bistable regime for \x > and a > noise-induced transitions 
between the metastable equilibria occur. 



analytically as {(x,/i) € IR 2 : x = 0} and {(x,fi) e M 2 : x = ±a//I} one has to use numerical 
techniques, such as numerical continuation, for more general systems. 

Interesting noise-induced dynamics occurs in the bistable regime for fi > 0. Fix any 
> 0, cr > and initial condition Xq. Then consider the first hitting times := inf{t > : 
x t = x^}. A standard result from probability [35] is that 

P(t± < oo) = 1 (4) 

i.e. no matter where we start, we will eventually visit both deterministically stable equilib- 
rium points with probability one. Although the result (jl]) is of importance from a theoretical 
viewpoint it is of very limited practical use. In particular, the time scale on which the stochas- 
tic switching between the potential minima occurs is of major interest. Suppose we start the 
process x t at Xq = x + . If a 3> 1 frequent switching occurs and we will quickly visit x~ while 
for < a ^ 1 switching is rare; see Figure |2J The theory of large deviations considers 
the first-exit time over the saddle point x* given by r + := inf{t > : Xq = x + , x t < x*} and 
shows that the mean first exit time is 

E[r+] = O LW^)-u^+)]/^\ ^a^Q. ( 5 ) 

The result ^ is also known as Arrhenius' law [B] and the rate 1/E[r + ] is called Eyring- 
Kramers rate [M1E3]; see also Section EO Furthermore observe that the potential difference 
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in (jSJ) is given by 



Hence the switching probability/rate also depends on the bifurcation parameter \x and in- 
creases when n — > + . In Figure [2] we show three time series for a fixed noise level with 
varying bifurcation parameter fi > over a fixed time interval t G [0, 1000]. It is clear that 
the dynamics in Figure EJ^b) with very frequent stochastic switching is different from rare 
switching events in Figure Ws) arid no switching events up to t = 1000 in Figure EJ^d). 
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Figure 2: (a) Bifurcation diagram with stochastic neighbourhoods of the metastable equi- 
librium points. The deterministic equilibrium curves are shown in black. The stochastic 
variance neighbourhoods B ± (h) defined by (I12j) are indicated in blue/green with confidence 
level h — 6. The overlap of B + (h) and B~(h) for fi = 0.5 is marked by a thick line (ma- 
genta). (b)-(d) Time series of (j2J) with noise level a = 0.5 for (b) \i = 0.5, (c) /i = 1.5 and 
(d) /i = 2.5. The dashed curves (red) indicate the metastable equilibria x . 



One possibility to capture the stochastic behaviour is to solve the forward Kolmogorov 
(or Fokker-Planck) equation [M] associated with <^ given by 

<9 d a 2 d 2 

—p(x, t) = —Q^(f(x] m(x, t)) + Yd^ p ( x > *) 

where p(x,t) = p(x,t\x ,t ) denotes the transition probability density of the stochastic 
process x t starting from xq at time to- However, solving flfjj) essentially solves the SDE (j2j) 
everywhere in phase space. It is clear that for higher- dimensional nonlinear problems - where 
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we are only interested in the local metastability of a equilibrium points or invariant sets - 
solving the PDE (JH]) may not be the best approach numerically. For small noise intensities 
- which are commonly assumed in applications - this is particularly unfortunate since the 
stochastic dynamics is very close to the zero noise limit a = on short time scales. 

Our approach tries to avoid these difficulties and aims at a natural extension of numerical 
continuation. We linearize fl2]) around the equilibrium points x which yields 

dX t = (/i - 3(a; ± ) 2 ) Xdt + adW t =: A(x ± ; fj)Xdt + adW t . (7) 

Observe that (J7J) is an Ornstein-Uhlenbeck (OU) process [38]. We will use the variance of 
the OU process to obtain neighbourhoods of x^ 1 within which sample paths of (j7|) stay with 
high probability. If the initial condition for (J7|) is deterministic then the variance of X t is 

Var(X 4 ) = a 2 [ u(t, s) 2 ds (8) 
Jo 

where u(t, s) is the fundamental solution [51] of the system 

u' = A{x k \n)u, u(t,t ) = l. (9) 
Defining V t := Var(X t ) direct differentiation of gives that V t satisfies the ODE 

V = 2A(x ± ;fi)V + a 2 . (10) 
Since \i > we have that A(x ± ; fi) = — 2/x < so that (ITU1) has a stable equilibrium point at 

Next, consider neighbourhoods of x ± given by the variance ffTTl) of the linearized process (see 
e.g. M 

B ± {h) := jxGM: \x - x ± \ < yjv(fi,(x)h = ^| (12) 

where Vv can obviously be interpreted as the standard deviation. The main idea of definition 
( JT2l) is that sample paths of (E]) stay with high probability inside B ± {h) if they are started 
at (or near) x ± . The parameter h scales the variance neighbourhood and can be used to 
control the probability to stay inside J B ± (/i) for a given time. Hence we can think of h as 
adjusting the confidence level of our metastable prediction (h = 1, one standard deviation; 
h = 2, two standard deviations; etc.). Figure EJ^a) shows B^ih) for three different values 
of /i = 0.5, 1.5,2.5 with fixed noise a = 0.5. This demonstrates that ( |T2|) can be used to 
approximate metastability properties for small noise intensities. Obviously all calculations 
for the SDE (J2]) can be carried out analytically. The open question is whether this approach 
can be used to construct a general and efficient numerical method. There are several problems 
that have to be considered: 

(PI) Generalize the construction of B ± (h) to arbitrary n-dimensional SDE. We summarize 
this well-known construction and the relevant results from probability theory in Section 

El 



6 



(P2) Find an efficient way to compute the covariance matrix of an OU-process during nu- 
merical continuation and/or for all points on a given equilibrium curve. The important 
step to solve this problem efficiently is to realize what information is already available 
from the deterministic continuation algorithm that can be used to compute the co- 
variance matrix. The main techniques from numerical analysis and control theory are 
summarized in Section HI 

(P3) Construct and efficiently compute a test function that detects overlaps of different 
neighbourhoods B ± {h). We suggest a test function based on the distance between 
ellipsoids. From computational geometry and optimization it is known that the dis- 
tance can be calculated by solving an optimization problem. The definition of the 
distance and all computational details are given in Section |5j 

Let us point out again that (P1)-(P3) are essentially all solved (or almost solved) as 
unconnected problems in various branches of numerical analysis, control theory, dynamical 
systems, optimization and probability. Our main contribution is to recognize the interplay 
between the different components which will provide a direct extension of deterministic 
continuation algorithms to metastable stochastic problems. 

3 Metastability and Linearization 

In this section we address the problem (PI) following Berglund and Gentz [IB] . Let x G M. n 
and consider the SDE 

dx t = f(x t ; fi)dt + aF(xt] n)dW t (13) 

where W t = {Wi,u Wb.tj • • • , ^k,t) T is standard k-dimensional Brownian motion, a > con- 
trols the noise level, /x G R is a parameter and / : M. n x R — > M. n and F : W 1 x M — > M. nxh are 
sufficiently smooth maps. Suppose the deterministic part of (fT3|) given by dx t = f(x t ;fi)dt 
has a hyperbolic stable equilibrium point x* = x*(fj,) for a given range of parameter values. 
Using a translation x = x — x* we get 

dx t = f(x* + x t ] fi)dt + o~F(x* + x t \ ii)dWf (14) 

Assuming that F(x*; /i) ^ the approximation of (|T4|) to lowest order via Taylor expansion 
is 

dX t = A(x*; n)X t dt + aF(x*; fx)dW t (15) 

where A(x',fi) = (D x f)(x; //) G M. nxn is the usual Jacobian matrix. Equation f[To"j) is an 
n-dimensional OU process. We assume that the initial condition xq is deterministic. The 
generalization of the variance © is the covariance matrix 

C t := Cov{X t ) = a 2 [ U{t, s)F{x*; fi)F{x*; fi) T U{t, s) T ds 
Jo 

where U(t, s) is the fundamental solution of U' = A(x*; jj)U . Differentiation shows that C t 
satisfies the ODE 

C = A{x*; fi)C + CA(x*; fi) T + a 2 F{x*; fi)F{x*; (if. (16) 
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Since x* is a hyperbolic stable equilibrium point, it follows [13j US] that the eigenvalues of 
the linear operator 

L(C) := A(x*; fi)C + CA(x*; /i) T 

are given by {2Aj}" =1 where Xj are the eigenvalues of A(x*; //) (and of A(x*; /i) T ). Therefore 
(Fl6|) has a stable equilibrium solution which is obtained by solving 

= A(x*; fi)C + CA(x*; + a 2 F(x*; fi)F{x*\ /i) T . (17) 

Observe that (TTTj) is a Lyapunov equation. It is well-known (see e.g. [51]) that the stability 
of x* implies the unique solvability of (TTT)) . For notational simplicity we shall not denote the 
solution of (TTTj) as C but simply write the symmetric covariance matrix as C or C(x*;fi). 
The main step of solving (I17p numerically at a given parameter value /i can be found in 
Section HI Then one can define a generalization of the variance neighbourhood from Section 
[2] as 

B(h) := {x E R n : (x - x*) T C~ l (x - x*) < h 2 } (18) 

where h is a parameter that can be interpreted as a probabilistic confidence level. A priori, 
the set (1181) may not be well-defined as C may not be invertible. It is well-known from 
control theory [109[ [99] that C is invertible if and only if the matrix 

Con(A, aF) := [aF oAF ■ ■ ■ aA n ~ x F\ G R nxnk (19) 

for A = A(x*; fx) and F = F(x*; /i) has maximal rank; this is sometimes concisely expressed 
as referring to the matrix pair (A(x*; fi), F(x*; /i)) as controllable |109[ [99] . From the con- 
trollability condition it follows that the invertibility of C is related to the structure of the 
noise encoded in F(x*]fi). In Section H] we discuss the case when C is not invertible. For 
now assume that C is invertible in which case the set 13(h) is immediately recognized as a 
solid ellipsoid with shape matrix 

Q := h 2 C, 13(h) = {x eW l : (x- x*) T Q-\x - x*) < 1} . 

It can be shown [TB] that stochastic sample paths stay in 13(h) near metastable equilibrium 
points with high probability. Similar results can also be found in the theory of large devi- 
ations [36]. It is quite lengthy to state the detailed asymptotic estimates depending on a, 
h and the eigenvalues of A(x*; fi). Since we are focusing here on a numerical algorithm we 
refer the reader to [16J for details. 

4 The Lyapunov Equation 

The next step is the numerical solution of the Lyapunov equation for a given metastable 
equilibrium x*(/i) as well as for an entire branch of equilibrium points obtained via con- 
tinuation 7 = {(x*(/i), fj,)} C M n x R. The algebraic equation (|T7|) is a uniquely solvable 
Lyapunov equation of the form 

AC + CA T + B = (20) 



8 



where we are going to use the shorthand notations A = A(x*; /i) and B := a 2 F(x*; fi)F(x*; jj) T 
from now on. Lyapunov equations have been studied in various branches of mathematics 
[3T] . Recall [54] that if one sorts the elements of C and B in vector form 

C T = (Cn,C 2 i, • • • ,Cl2, • • ") T and b T = (611,621, ••• ,&12,"--) T 

then (1201) can be rewritten as a standard linear system 

[I®A + A®I]c=-b (21) 

where ® denotes the Kronecker product [12] • The problem of efficient numerical solution of 
(120]) or (1211 (and of several generalizations) has attracted considerable attention in numerical 
analysis and control theory [37]. For our situation several new aspects arise since we want 
to solve f )20|) along an entire equilibrium branch 7: 

1. All standard numerical continuation algorithms require an approximation of the nx(ri+ 
1) Jacobian matrix (Df X)l Af)(x*(/j l x), fix) to compute a point (x*(/i 2 ), /i 2 ) G 7 starting 
from (x*(/ii), fi±) G 7. Therefore, the matrix A = (D x f)(x*(^i); fix) is available at each 
continuation step. Furthermore, computing the matrix B requires at most one matrix 
multiplication at a given point (x*(/ii), /ii); for purely additive noise B can even be 
precomputed for all equilibrium points. 

2. Solving (120]) at (x* (fix) , fix) gives a matrix C(x*(/ii); /ii). If |/ii — /i 2 | is small then 
C(x*(/ii); /ii) is already an excellent initial guess to find C(x*(/i 2 ); fi 2 )\ Hence, except 
for the first point on the equilibrium curve, we always have an initial guess available 
for iterative methods. 

The observations suggest that computing the covariance C should be relatively easy. We 
decided to focus on three different approaches which we briefly review here. Due to a good 
initial guess, the most natural choice are iterative methods. Consider the reformulation ff2T]) 
and define A := [/ <g) A + A <g> I] . Then the standard Gauss-Seidel iteration [112] is given by 

c (fc+i) = c (fc) _ M^ s (Ac {k) + 6) (22) 

where Mqs is the matrix obtained from A by setting all entries above the diagonal (Aij, 
j > i) to zero. The iteration is terminated when ||c^ fe+1 ^ — c^|| < tol where tol is a given 
tolerance. Other possibilities for iterative methods include the Jacobi method and successive 
overrelaxation (SOR) methods [112] . For large sparse Lyapunov equations several special 
methods have been suggested including alternating-direction-implicit (ADI) by Wachspress 
[114] and special SOR methods by Starke [111] . We shall not consider the special methods 
here although they should definitely be relevant for large scale bifurcation problems [1UU] . 

Another well-known method for the iterative solution of (120 p is Smith's algorithm [106] . 
The first step is to fix a scalar q > and consider the matrices 

K := 2q(ql - AY l B(qI - A T Y 1 , 
G := (ql - A)-\qI + A). 
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Direct matrix multiplication shows that (|20|) is equivalent to solving 



C = K + GCG T . (23) 

The iteration of (|23|) converges linearly. Smith observed that with initial guess = K 
the iteration 

C (k+i) = C (k) + G 2<° C (k) r G 2*y ; fe = o, 1,2, ... (24) 

obtained by squaring G at each step converges quadratically. The algorithm is terminated 
when HC^ 1 ) — C^|| < tol. Using ADI theory the optimal q > can be found and the 
error has been calculated |114[ [37]; we will simply fix q = 0.1 which is the classical choice by 
Smith |107j . Observe that Smith's algorithm does not use an initial guess. 

There are also several direct (non-iterative) algorithms available. The most important 
techniques were suggested in the 1970s [HI E2J HI] and have become standard methods for 
the numerical solution of (12"U|) . The Bartels-Stewart algorithm [9] requires to compute the 
real Schur decomposition of A given by U T AU = R where U G R nxn is orthogonal and 
R g M. nxn is upper quasi-triangular (i.e. diagonal with possible 2x2 blocks on the diagonal 
corresponding to complex eigenvalues). Then (1201) can be transformed to 

RY + YR T = W (25) 

where W = —U T BU and Y = U T CU. The right-hand side W can be obtained by solving 
UW = BU for W . Solving (^o§ requires the solution of an upper quasi-triangular system 
which is straightforward. Then one can solve CU = UY for C to get the final result. The 
Bartels-Stewart algorithm can also be helpful for our problem as it can be used to solve the 
problem at the first continuation point and it applies when Gauss-Seidel iteration fails as a 
"fall-back" strategy. 

In Section [7J we are going to compare the performance of the Bartels-Stewart algorithm, 
Smith's method and Gauss-Seidel iteration for a practical numerical continuation problem. 
Once we have the covariance matrix C it is important to check whether C _1 exists so that 
()12p is a well-defined ellipsoid. We are going to illustrate why such a test is important. In 
control theory [69] it is well-known how to define ellipsoids in a degenerate case when the 
shape matrix Q = h 2 C is only positive semidefinite (see also Section ED 

B(h) := {xeR n : v T x < v T x* + (v T Qv) 1/2 G W 1 } . (26) 

Now consider the example 

A = A(x*; fi) := ^ ~ 2 ^ ) , B = a 2 F(x*; f,)F(x*; := ( ^ ° 

which corresponds to a stable hyperbolic equilibrium in M 2 with additive noise on the first 
component only. Even without solving for C we can compute the matrix ffT9l) 



Con(A(x*; ji), aF(x*; //)) 



a -2d 
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so that C is not invertible since rank(Con(A(x*; /i), aF(x*; //))) = 1. Indeed, we easily find 
that solving the Lyapunov equation gives 

Assuming for simplicity that x* = (0, 0) T we get that the set 13(h) defined in ( 126]) is given 
by 

{/ v 2 2U2 ] 
(x u x 2 ) G M 2 : v x xx + v 2 x 2 < V(vi, v 2 ) G M 2 ^ 

= {(xi,x 2 ) G M 2 : < crft/2, x 2 = 0} . 

The ellipsoid 13(h) is a degenerate interval which reflects that the degenerate noise terms only 
act on the xi-coordinate. Detecting such a degenerate (or near-degenerate) noise is clearly 
important in applications as this identifies directions along which metastable escapes are 
unlikely. A simple test for this degeneracy is to compute the singular value decomposition 
(SVD) 02] of C. 



5 Ellipsoids and the Testfunction 

Suppose we have two covariance matrices Ci j2 for given set of parameter values at (/i, x\ 2 (a0) = 
(/i, x\ 2 ) and we are interested in detecting the distance between the associated ellipsoids as 
large/small distances are expected to correspond to long/short travel times of sample paths. 
Denote the shape matrices of the ellipsoids by Qi j2 = h 2 C\^ and a general ellipsoid by 

£ = £(x*, Q) = {xeR n :(x- x*) T Q~ x (x - x*) < 1}. 

The idea of considering covariance ellipsoids and their overlaps is not new. For example, 
the idea is used in satellite tracking for collision avoidance [2]. In computational geometry 
and robotics one often considers the minimum-volume enclosing ellipsoid of an object, also 
called the Lowner-John ellipsoid [94] . In statistics an ellipsoidal distance defined via the 
covariance matrix, the so-called Mahalanobis distance |76] , is often used [96] • Various method 
have been proposed to detect ellipsoid overlaps ranging from Grobner bases [2T], analytical 
representation formulas |3], reformulation as an eigenvalue problem [M], local approximation 
by balls |74) to polyhedral approximations jlOj [18]. Here we will adapt an idea based on 
calculating the distance between ellipsoid by solving an optimization problem which has 
several advantages to be discussed below. The support function of an ellipsoid is [69] 

(Js(v) := supt> T x = v T x* + (t^Qf) 1 / 2 . 

The Hahn-Banach Theorem [7] gives that an ellipsoid with a positive semi-definite shape 
matrix can be defined as 

S := {x G R n : v T x < v T x* + (v T Qv) 1/2 \/v G R n } . (27) 
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A measure of the distance between two ellipsoids |70j is given by 
6 = 8(£(x* 1 ,Q 1 ),£(xl,Q 2 )) = max {-cr £l (-v) - ae 2 (v)) 

{v T x{ - {y T Qiv) l/2 - v T x* 2 - {v T Q 2 vf' 2 ) . (28) 



max i V X-i 



The distance 5 and its definition have several advantages for detecting metastability. The 
definition also applies immediately if the matrices Qi j2 are degenerate. For example, if we 
are interested in the distance of a covariance ellipsoid £(x*,Qi) to an unstable equilibrium 
point x* 2 (e.g. the saddle point in Section [2]) we can just set Q 2 = and still consider the 
distance 5. We can even replace the ellipsoid with a more general convex set % if the support 
function <Jn( v ) is easy to calculate. The main advantage is that S is also a test function since 

• 5(£ (x\, Qi), £(x 2 ', Q 2 )) > if the two ellipsoids are disjoint, 

• 5 (£(#*, Qi), £(x 2 , Q 2 )) = if the ellipsoids touch at a point, and 

• 8(£ (xi, Qi), £{x\, Q2)) < if the ellipsoids intersect. 



Therefore the distance (1281) is a test (or bifurcation) function if we want to check how 
likely metastable transitions occur in our SDE (fl3|) . Note that the precise number of noise- 
induced transitions cannot be inferred from 8 as we have not made any assumptions about 
global dynamics; but see Section |9j Observe that ( 128]) is a classical nonlinear optimization 
(or nonlinear programming) problem. In standard minimization form with a differentiable 
constraint it can be written as 

min {—v T x* 1 + (v T Qiv) 1/2 + v T x* 2 + (v T Q 2 v) 1/2 ) =: min G(v), , , 

subject to = ||-u|| 2 — 1 =: g(v). 

and we obtain a solution to ( 128]) by the negative solution value of ( 129]) . Many efficient 
algorithms for the solution of ( 129]) are available [83]. In particular, many iterative schemes 
are known among which sequential quadratic programming (SQP) [52[ [90] has turned out 
to be among the most powerful techniques. Here we simply use this approach which solves 
a quadratic programming problem at iteration step k given by 

min [\w T H k w + VG(v k ) T w) , , , 

subject to Vg{v k ) T w + g(v k ) = 0, 



where H k = (y 2 L)(v k ,Uh) is the Hessian of the Lagrangian L(v,u) := G{v) — u T g{v) and 
u k G R is an approximation of the Lagrange multiplier. If the solution of ( 130]) at step k 
is denoted by w k then the main iteration step is v k+ i = v k + a k w k for a given step length 
a k > 0. It is important to note that VG(v k ), Vg(v k ) and H k can be supplied in explicit 
form 

vgk) = -xt + ^(^Qi^)- 1/2 (gi + gD^ + ^ + ^(^2^)~ 1/2 (g 2 + QD^ ; 

Vg{v k ) = 2x k , 



{(#*)y}?j=i = -i(vlQxv k )- 3/2 [(Qi + Ql)v k U(Qi + QDv^ + -{vlQ x v k )[Q x + Ql 



-^Q 2 v k )- 3/2 [(Q 2 + Ql)v k ]i[(Q 2 + QDv^ + -{vlQ 2 v k )[Q 2 + Q£] - 
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which avoids the computation of finite difference approximations during the optimization 
iteration. We use a standard quasi- Newton line-search method to solve ( I30p . The iterative 
algorithm stops when the solution, solution values and constraints are below a given toler- 
ance. The details of this part of the algorithm will not be discussed here and details can be 
found in [831 [78]. 

As for the Gauss-Seidel method, it is very important to point out that the iterative 
solution of fl2B"j) can be used efficiently during continuation. Given a fixed point x*(fi\) at 
parameter values \i\ we obtain a solution v(ni) to (1281 by solving ( 1251) . For an equilibrium 
point continuation step from Hi to /i 2 we have that |/ix — fi^l is small so that v (/ii) can be 
used as a very good initial guess for the optimization problem to be solved with parameter 
values /i2. 

6 Algorithm Summary 

In this section, we summarize the main steps of our algorithm which augments deterministic 
numerical continuation. Consider the SDE 

dx t = f(x t ; n)dt + aF(x t ; ^)dW t , for x t G R n and /i G R. (31) 

We assume that a stable equilibrium x*(/j i0 ) for the deterministic part of f l3T|) is given (or it 
can be found e.g. using Newton's method [112J ) so that f(x*(no); /i) = 0. Then define 

A(x*(ji )\fi ) = (D x f)(x*(fj, ); fi ). 

Using the Bartels-Stewart algorithm (see Section H|) we solve 

= A(x*(fi ); Mo)C + CA(x*{hq)- /i ) T + a 2 F(x*(fi ); /i )F(x*(/i ); fi ) T 

for the covariance matrix C = C(x*(/j,q); /j,q). This completes the initialization step. The 
main iterative step of the algorithm is as follows: 

(Al) Choose a step length [3 k and set //& = Hk-i + 0k- Solve the continuation problem 
f(x;Hk) = for the new equilibrium x*(nk) with starting point x*(fik-i) (see e.g. 

hid. 

(A2) Consider the Lyapunov equation 

= A(x*(fi k ); n k )C + CA(x*(fi k ); /i fc ) T + a 2 F(x*(fi k ); fi k )F(x*(/i k ); fi k ) T 

and solve it for C, preferably using an iterative algorithm with initial guess C(x*(fJ,] t -i)', A*fe-i)- 
This yields the new covariance matrix C(x*(fi k ); fi k ). Define the shape matrix Q(x*(fi k ); fi k ) := 
h 2 C(x*(n k ); n k ) for a given confidence level h. 

(A3) Given two shape matrices Qx.fe := Q{ x \{^k)'i ^k) an d Qi,k '■= Qix^ifik); fi k ) for dif- 
ferent stable equilibria, compute the distance S = S(fi k , h) between the two ellipsoids 
£i(x*(fik);Qi t k) and ^(a^O^fc); Qi,k) solving the optimization problem ([28]) using an 
iterative method such as SQP (see Section [5]) with initial conditions obtained from the 
iteration step k — 1 . 
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As an output we get the following parameterized families 

• equilibrium points from numerical continuation, 

• ellipsoids {£ (x*(/jl), C(x*; from solving Lyapunov equations, and 

• mutual distances {S(/j,, h)}^ from solving a nonlinear programming problems. 

The ellipsoids {£(x*(fi), C(x*; provide locally rigorous estimates for metastability 

[16]. The distance 6(f/,k, /i) between two ellipsoids an d Q2,fc gives an indicator for global 
transitions occurring from Qi.fc to Q2,fc or vice versa, based on the assumption that larger 
distances correspond to lower switching probabilities. Section [7J shows that using h) 
works nicely in practice. Nevertheless, it may be desirable to obtain rigorous estimates if 
global assumptions are made; Section [9] augments the algorithm in this direction for a special 
case. 

Furthermore, it looks intuitive to consider higher-order moments of the fully nonlinear 
stochastic process described by the SDE ( {TBI . However, the ODEs for higher-order moments 
usually do not form a closed system [38j such as ( jl~5i) . Observe carefully that if a set of 
moment equations forms a finite-dimensional closed system (or can approximately be closed 
[108J) then a modified version of steps (Al)-(A2) should carry over to this situation since 
equilibria for the moment ODEs satisfy an algebraic equation which can again be solved 
iteratively with a good initial guess from the previous continuation step. 

We note that the algorithmic steps (A2)-(A3) can be used as a post-processing tool for 
an existing numerical continuation curve (x*(fi); fx). This is not as efficient as combining 
(Al)-(A3) as it requires re-building the matrices A(x*;fi). In summary, we have obtained 
local approximate information about a system of stochastic differential equations using a 
completely deterministic continuation algorithm. The additional computations required to 
obtain this information are easy to implement in a classical continuation algorithm and/or 
bifurcation software package. The iterative solution procedures for the Lyapunov equation 
and the ellipsoid distances are expected to make the algorithm computationally very efficient. 

7 Neural Competition and Bistability 

In this section we are going to test our algorithm for the situation where the deterministic 
dynamical system has two stable coexisting equilibrium points ("bistability"). The differ- 
ential equations we are going to study are based on ODEs modelling a two-cell inhibitory 
neural network [2U [23] . The goal is to describe competition between two neural populations. 
For example, such a situation can occur due to ambiguous external stimuli [73] inducing a 
bistable behaviour in the neuronal system. A typical example is binocular rivalry [31] where 
switching between different visual perceptions occurs. This situation can be modelled [104J 
by considering the (fast) spatially averaged firing rates X\^ of two neural populations and 
two associated (slow) time fatigue accumulation variables y\ %2 - The resulting ODEs are 




-xi +S(I C - [3x2 - gyi) 

-x 2 +S(I C - 0xi - gy 2 ) 

- Vi), 
e(x 2 -y 2 ), 



(32) 
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where I c is the main bifurcation parameter and the sigmoid-shaped gain function S : R — > R 
is often chosen [21] in numerical simulations and continuation calculations as 



l + exp(-r(u-0))' 

We adopt this choice and also fix the parameters 

= 1.1, # = 0.5, r = 10, 9 = 0.2 (33) 

so that our calculations are a direct extension of numerical continuation in [23] • The param- 
eter < e < 1 describes the time scale separation between the fast and slow variables. We 
are only going to consider fl32|) in the singular limit e = of perfect time scale separation. 
The equations 

x\ = -xi +S(I C - /3x 2 - gyi) =: fi(x), ^ 
x' 2 = -x 2 + S(I c - /3x 1 -gy 2 ) f 2 (x), 

are also called the fast subsystem of (1321) where y\ i2 are regarded as parameters. For an in- 
troduction to the theory of fast-slow systems and singular limits see [561182]; an example how 
fast subsystem bifurcation analysis can form a building block of bifurcation analysis for the 
case e > can be found in [48, 49J. Since (I34j) is a model for the activity of (finite) neuronal 
populations there are various natural stochastic effects such as channel noise |35J , input noise 
[113] , neuronal background noise |60j and external noise in experiments/observations [57J. 
Therefore it is reasonable to extend (EH) to the SDE 



(i Xl ) = ( ~ Xl + S 3 ~ " m \ ) dt + a 2 F{x)dW t (35) 
\dx 2 J \ -x 2 + S{I C - (5xi - gy 2 ) J 

where F : R 2 — > R 2x2 . Furthermore, we fix the slow variables to 

y x = 0.7 and y 2 = 0.75 (36) 

which introduces a slight asymmetry into the system. Both slow variables also lie within 
plausible ranges as considered in [24] . 

Figure [3] shows a continuation calculation for the neuronal competition model fl35l) with 
parameter values fl3"31 and f[3"6"l) . The additive noise terms are given by 

a 2 F(x*)F(x*) T = a 2 ( ^ °f ) for a = 0.3. (37) 

The deterministic equilibrium continuation has been carried out using the Moore-Penrose 
algorithm [7T], |44] with fixed continuation step size 0.001. For the the computation of the 
covariance ellipsoids and the distance between them we refer to the summary of our algorithm 
in Section El Figures E^a) and (d) visualize the ellipsoids and Figure MJo) shows the distance 
between the ellipsoids defined by (|2"8"]) . We find two regions where the distance is negative 
and overlaps between ellipsoids occur. Hence we expect that equilibrium points in the 
parameter regions with overlaps are only weakly metastable and relatively frequent noise- 
induced switching between different neuronal activity patterns occurs. This conjecture is 
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Figure 3: Continuation results for (1351) with parameter values (I33p and (I3"6"|) . The noise terms 
are given in (1371) . (a) The thick curves (red and blue) show stable equilibrium point branches 
continued in the main bifurcation parameter I c . There are two saddle node bifurcations on 
the isolated bifurcation curve (isola) from a stable node (blue, thick curve) to a saddle (green, 
thin curve). We also show some of the two-dimensional ellipsoids £\^ calculated in (21,0:2) 
phase space for fixed parameter values I c and embedded in (xi, X2, /^-coordinates; here we 
use h — 1 as the confidence parameter in the definition of the covariance matrix, (b) Calcu- 
lation of the distance S = S(£\, £2) between the ellipsoids in the bistable regime as defined in 
( 128|) . (c) Mean number of passages T p of a trajectory between the two stable equilibria over a 
time interval [0, 1000]; averaged results over 100 sample paths obtained via direct numerical 
integration of the SDE are shown (grey, dots indicate grid in I c ). (d) Projection of (a) onto 
(J c , X\) where the X\ maxima and minima of the ellipsoids have been connected to form tubes 
around the stable equilibrium branches, (e)-(f) Direct numerical SDE simulation for I c = 0.7 
(as indicated by the arrow from (d)). The colored dots correspond to the equilibria. 



confirmed in Figure [3](c) where the mean number of noise-induced passages T p between two 
stable equilibrium points pi^ is shown during a fixed time interval; more precisely, consider 
a trajectory j(t), fix some small p > and define 

T^ 2 (l) := h) :t 1 <t 2 <T max ,\\ 1 (t 1 )-p 1 \\ < p, t 2 = inf{t : t > h, \\j(t) - p 2 \\ < p}} 

T^\l) := #{{h,h) :t l <t 2 < T max , || 7 (ti) -p 2 || < p,t 2 = mf{t : t > t u || 7 (t) - Pl \\ < p}} 

which just count the number of times a trajectory starting from a small ball near pi reaches as 
small ball near p 2 and vice versa. Then we can average the results over different realizations 
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of the noise (i.e. over different paths 7) 



(3? 



For Figure E^c) the parameters p = 0.05 and T max = 1000 have been used and the expected 
value in (|38p has been computed over 100 sample paths. Note carefully that distance function 
in Figure Mb) predicts the qualitative shape of the passage time distribution T p very nicely. 
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Figure 4: Computation of the covariance matrices along the entire two stable equilibrium 
continuation curves (see Figure |3]) for given tolerances (logarithmic abscissa). The colors 
(red/blue) indicate the equilibrium branch as in Figure [3j (a) Total computation time (in 
seconds, linear interpolation of times is shown) for different tolerances of \\C^ — C( fc_1 )|| < 
tol. Gauss-Seidel iteration (crosses), Smith iteration (stars) and a direct solution via the 
Bartels-Stuart algorithm (as implemented in [78]) are compared, (b) Comparison of the 
total number of iteration steps for the Gauss-Seidel (crosses) and Smith (stars) algorithms 
are compared. 

We shall not investigate the dynamical implications from our method here but focus on 
the performance of the algorithm. As a starting point we use the two continuation curves 
of stable equilibrium points shown in Figure [3j For each curve we calculate the covariance 
matrix by solving the Lyapunov equation for each point on the continuation curve using 
Gauss-Seidel and Smith iterations as well as the Bartels-Stewart algorithm. For the Gauss- 
Seidel algorithm we use as the starting point of the iteration the covariance matrix from the 
previous point on the equilibrium curve. Figure H] shows the computation time as well as 
the average number of iteration steps along the equilibrium curve for different tolerances of 
the iteration termination condition 

We see that for relatively low tolerances between 10 -2 and 10~ 7 the iterative solution using 
the Gauss-Seidel method seems to perform best. This is not surprising since it is the only 
method that uses the previous point on the curve of equilibria which is expected to be 
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an excellent initial guess. For higher tolerances and high-precision computation Smith's 
algorithm as well as the exact Bartels-Stewart method seem to be preferable. Since Smith's 
algorithm always converges quadratically this is again expected in comparison to Gauss- 
Seidel. Using SOR or ADI iterative techniques or considering larger systems could potentially 
even further increase the advantage of iterative methods that use an initial guess from the 
previous point on an equilibrium curve; see also Section HI Another important conclusion 
from the calculations in Figure H] is that even though the two equilibrium curves have O(10 3 ) 
points each, the calculation took only a few seconds. Therefore the computation of all 
covariance matrices of equilibrium curves is expected to very fast on standard single-machine 
computer hardware for most small to medium-size ODE systems. 
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Figure 5: Computation of the distance between the two covariance ellipsoids of two stable 
equilibrium continuation curves (see Figure |3]) for given termination tolerances (logarithmic 
abscissa) of the SQP optimization algorithm. The distance is calculated at every tenth 
point equilibrium point; 216 distances have been computed in total, (a) Average number of 
function evaluations of the objective function over the 216 points required during the SQP 
algorithm, (b) Average number of iteration steps, (c) Total time in seconds to process 216 
points. 

Figure |5] shows an overview of the computational cost to obtain the distances shown in 
Figure [3]^ c) between ellipsoids using SQP as implemented in [7S]. The distance has been 
computed for 216 covariance ellipsoids sequentially along the equilibrium point curves. The 
initial conditions were obtained from the result of the previous optimization problem. The 
main result of Figure [5] is that the distance calculation can be carried out quickly and re- 
quires very few iterations steps and function evaluations. This means that we can evaluate 
the testfunction for overlapping ellipsoids efficiently using optimization. However, we do 
not claim that the algorithm we used here is optimal in any way. It is possible that other 
optimization techniques of methods to estimate distances between ellipsoids outperform the 
SQP approach we used here. However, from a practical point of view the results we obtain 
show that the computational time is certainly not prohibitive to process entire equilibrium 
bifurcation branches. 
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8 A Predator-Prey System 



In the previous section, we have focused on the algorithmic cost of our algorithm and the 
distance calculation between ellipsoids. In this section we are going to consider an example 
with a complicated noise term and focus on the value of our method for applications. The 
classical Rosenzweig-MacArthur [98] model for the interaction of predators Y and prey X is 
given by 



where (x, y) represents the population densities of (X, Y), 7 relates to the carrying capacity 
of the prey, (3 is a conversion factor and m a parameter describing mortality of the predator. 
The model ( 139|) can be derived as a large-system size limit for the individual interactions 
between X and Y . Finite-size effects of the population can be included into a stochastic 
fluctuation term. Using a Kramers-Moyal (or system-size) expansion one finds [95, 38J 



where W t = (W} ,W} ) T is standard Brownian motion, the matrix- valued function C is 
given by 



and a = 0(1 /N) where N is the population size. Therefore o — > corresponds to the 
limiting case of an infinite population which recovers the deterministic limit ( |39l) . Observe 
that the noise terms in ( HOI) are multiplicative and exhibit correlations between the two pop- 
ulation densities. Therefore it is not immediately clear how a bifurcation diagram of ( 1391) is 
altered once the (more realistic) finite-system size is considered. 

We focus on deterministic Hopf bifurcations in the model which have received the most 
attention in ecological predator-prey models [62J . Figure El shows an equilibrium contin- 
uation in 7 where increasing 7 can be interpreted as increasing the carrying capacity for 
the prey. Observe that a stable focus undergoes a Hopf bifurcation. Classical deterministic 
ecological theory [97] argues that increasing the carrying capacity corresponds to enrichment 
and that the periodic solutions born in the Hopf bifurcation can move the system close to 
the coordinate axes 

{(2, y) e R 2 : x > and y = 0} and {(2, y) G K 2 : y > and x = 0} 

which delimit the positive quadrant. Once the system reaches any of the two axes it is easy 
to see that this corresponds to extinction of a species leading to a "paradox of enrichment" . 
This "paradox" is a highly debated topic in ecology and many different ways of resolving it 
have been suggested, see for instance [U [551 [391 ETJ- 




(39) 



y' = -my 




(40) 
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Figure 6: Continuation for ( HUj) in 7 with m = 1 and (3 = 3 fixed, (a) and (b) are two 
viewpoints fo the same bifurcation diagram. The thick blue curve is computed via equilibrium 
continuation. The thin blue ellipses are computed with the algorithm from Section [6]for h = 1 
and a = 0.01. The Hopf bifurcation (H) at 7 = 2 is marked with a green dot. The red planes 
delimit the positive quadrant in (x, y)-phase space. The black dots on the ellipse at 7 = 1.9 
are interpolation points for the ellipse outside the positive quadrant. 

However, from our computation the ellipsoids suggest a very simple solution. The 
predator-prey system before a Hopf bifurcation can easily reach the well, even for 

small noise which corresponds to a large (but finite!) population size. Close to the bifur- 
cation point the ellipsoids increase in size which is precisely the well-known slowing down 
effect exploited in the theory of critical transitions [66| Wf\ 1101] . If the carrying capacity 
in an ecosystem only increases slowly, which is reasonable to assume, then we expect that 
stochastic effects drive the system to extinction before the "paradox of enrichment" Hopf 
mechanism becomes relevant i.e. one would not see regular oscillations before extinction. 
Furthermore, the deterministic periodic solution occurring due to enrichment could actually 
have a stabilizing effect as the stochastic effects are small for a strongly attracting determin- 
istic periodic orbit far from bifurcation. Indeed, the idea of stabilization of enrichment has 
been considered previously [7H1 EB] ■ 

9 A Special Case - Kramer's Law 

So far, all computations only required local assumptions on the SDE f lT3|) regarding existence 
of a deterministic equilibrium and suitable smoothness. The covariance neigbhourhood 13(h) 
provides a rigorous local control of the dynamics. The global distance 5 is a precisely 
computable, but probabilistically heuristic, measure to gain insight into global transition 
dynamics; without global assumptions on the dynamics this seems to be the best we can 
hope for. However, one may ask what happens if we have additional information on the 
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global dynamics. Consider the SDE 



dx t = -VU^x t )dt + adW t , x G R n , U^.R n -> R, W t G W 1 , aeR (41) 

where the deterministic part is a gradient system with a potential parameterized by 
H G R. Critical points of £/ M correspond to equilibria for the deterministic dynamics. Fix 
some ft G 1R and suppose has precisely two local minima x* and y*, corresponding to 
stable equilibria, and one saddle point z*. Define 

Ttf( y *) := inf{t > : x t G N(y*)}, x = x*, 

for a suitable neighbourhood Af(y*) of y*. Under the assumption that the saddle point z* 
has a single unstable eigendirection with eigenvalue X(z*]fi) > 0, the Eyring-Kramers law 
EH [63 1 states that 



(42) 

where A(x*;/x) = D 2 U^ G M™ xri is the Hessian of and J\f a 2/ 2 (y*) a ball of radius a 2 / 2 
around y*. The precise formula (j4"2"]l is due to Bovier et al. [19]; see also [531 EH] for reviews 
and generalizations of Kramers' law. The probability of switching due to noise from x* to 
y* is given to leading-order by (|42|) and interchanging the roles of x* and y* provides the 
noise- induced switching estimates for the transition from y* to x*. To compute (|42p we can 
follow an analogous strategy as for the more general case discussed so far. The equilibria 
x*, y* and z* as well as the associated linearizations A(-; fi) can be computed efficiently via 
numerical continuation for a curve parameterized by fi G R, the function V is available by 
assumption and it remains to compute X(z*; fi), det(A(z*; ft)) and det(A(x*; fi)). To compute 
the determinants we can simply use the LU decomposition [42J which also works well for 
large sparse systems. However, calculating the leading eigenvalue is bound to be costly if 
we look to compute all eigenvalues and then extract the leading one. Suppose we are given 
the results A(z*; fik-i), X(z*;fik-i) and the associated eigenvector v(z*;fik-i) from at the 
last continuation step then we can again use an iterative method to compute X(z*; fik). For 
example, setting v (z*; fik) = v(z*; fik-i) an d Xq(z*; fi k ) = X(z*; fik-i) then Rayleigh quotient 
iteration [SB] is given by 



/ * \ {MA Vk) - \j{z*\ fi h )I) l Vj{z*\ Mfc) . ( . 

"' + ' (2;w)= ||(^.; W )-A i (,-; W )/)-^-; W )ll' fa ' = - 1 - 2 -- < 43 > 

and the eigenvalue for the j-th iteration step is 

w . n _ Vj( z *'> Vk) T A(z*; fM k )vj(z*; fi k ) (AA , 
Aj[Z ^ k) ~ v^fikYv^fik) ■ lMj 

It is well-known that for a symmetric matrix A(z*; fij) the iteration (14"3"1)-(14"4"1) converges cubi- 
cally to the leading eigenvalue and eigenvector [861 HO]; in particular, Xj(z*; fik) — > X(z*; fik). 
Since A(z*;fij) = D 2 U N is derived from a, sufficiently smooth, potential U H the matrix 
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at each continuation step is symmetric and the fast convergence results for Rayleigh itera- 
tion apply; note that this may not be the case for open sets of "bad" starting conditions 
if the matrix is not symmetric [Uj. In any case, evaluating the remaining terms in (1421) is 
straighforward so that we can calculate mean-first passage times between equilibria in gra- 
dient systems quickly, with high accuracy, and rigorous error estimates by using numerical 
continuation. 
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